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We study the Bloch dynamics of a quasi one-dimensional Bose-Einstein condensate of cold atoms 
in a tilted optical lattice modeled by a Hamiltonian of Bose-Hubbard type: The corresponding 
mean-field system described by a discrete nonlinear Schrodinger equation can exhibit dynamical 
(or modulation) instability due to chaotic dynamics and equipartition over the quasimomentum 
modes. It is shown that these phenomena are related to Bogoliubov's depletion of the Bose- 
Einstein condensate and a decoherence of the condensate in the many-particle description. Three 
types of dynamics are distinguished: (i) decaying oscillations in the region of dynamical instability, 
and (ii) persisting Bloch oscillations or (iii) periodic decay and revivals in the region of stability. 

PACS numbers: 



I. INTRODUCTION 

Recently, considerable attention has been paid to the 
dynamics of cold atoms and Bose-Einstein condensates 
(BECs) loaded into an optical lattice with a static, e.g., 
gravitational, force (see [l], Q and [!, 0, H, H, 0, H, H, 
Hoi . [Til . [TH ) . In the limit of vanishing particle interaction 
the system shows Bloch oscillations (BO) (for a review 
see, e.g., [13 )■ However, it is known that the interaction 
between the atoms leads to modifications and possibly 
even to a breakdown of these Bloch oscillations. 

The most prominent theoretical approach to investi- 
gate the dynamics is the reduction of the many-particle 
system to a mean-field description via the (nonlinear) 
Gross-Pitaevskii equation. The description can be fur- 
ther simplified by discretizing it in terms of Wannier func- 
tions localized on the lattice sites. In the many-particle 
description this yields a Bose-Hubbard model and ac- 
cordingly the mean-field dynamics is described by the 
discrete nonlinear Schrodinger equation (DNLSE) (see, 
e.g., [lij]). It should be pointed out that the mean-field 
approximation is formally equivalent to a classical limit 
for single particle quantum mechanics. Therefore, it is 
often denoted as '(pseudo) classical', although it still de- 
scribes a quantum system. This becomes evident in the 
limit of vanishing interaction, where it reduces to a single 
particle Schrodinger equation. Nevertheless, the formal 
similarity of the many-particle to mean-field transition 
and the quantum classical correspondence allows for the 
application of semiclassical methods as well as the in- 
vestigation of topics such as quantum chaos within the 
framework of ultracold atoms [H QJ, [T3, M, 03 M, HI • 

The mean-field system of interest in the present study 
shows a rich structure of mixed and chaotic behavior and 
our main aim will be to identify the counterparts of some 
pronounced features in the corresponding many-particle 
system. For this purpose we compare the DNLSE dy- 
namics with the underlying microscopic many-particle 



system, described by the ID bosonic Hubbard model with 
the Hamiltonian 
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where aj and di are bosonic creation and annihilation 
operators for the Ith. lattice site, and hi = a\ai are the 
associated number operators. The hopping energy and 
the particle interaction are denoted by J and W, respec- 
tively, and ei = dFl is the on-site energy where d is the 
lattice period and F the magnitude of a static force. Ex- 
perimentally, all these quantities can be controlled sepa- 
rately. However, the validity of the Bose-Hubbard model 
is based on certain assumptions such as the single band 
tight binding approximation, which may break down in 
some parameter regions. Yet, for the following studies 
only the ratios of the quantities are of interest for the 
qualitative behavior of the system which offers an ad- 
ditional freedom to maintain the validity of the Bose- 
Hubbard Hamiltonian. Therefore, we believe that the 
observations reported in the following should be experi- 
mentally accessible. 

The Hamiltonian commutes with N = hi and 
therefore the total number N of particles is conserved. 
Here we consider the macroscopic limit N — > oo, W — > 
with constant g = WN/L where the number of sites, L, 
is kept finite. In this limit the mean-field approximation 
can be applied, which usually is formulated as replacing 
the bosonic operators ai, aj by complex numbers ai, ct[, 
the components of an effective single particle wave func- 
tion which appear as 'classical' canonical variables. The 
resulting mean-field Hamiltonian function is given by 

=E e <N 2 -^E ( a ^ ai + c-c)+f£N 4 ' (2) 
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up to a term proportional to | cz; | 2 , which is an inte- 
gral of motion. The DNLSE can be formulated via the 



2 



canonical equations of motion 

ihai = dH/dai , iHal = —dH/dai. (3) 

One aim of the present paper is to investigate the valid- 
ity of this mean-field approach for a system with finite 
particle number in the different parameter regions. 

For nonvanishing particle interaction g ^ the mean- 
field Bloch oscillations are damped [1, Q. In this con- 
text the main phenomena are the dynamical instabil- 
ity which is also known as modulation instability (see, 
e.g., [22l I23I ]) and the equipartition over the quasimo- 
mentum modes (also denoted as thermalization) due to 
the onset of classical chaos in the DNLSE. On the other 
hand, within the many-particle approach the main phe- 
nomenon induced by the interaction is a decay of Bloch 
oscillations due to decoherence @, HH . The present anal- 
ysis shows a direct relation between these classical and 
quantum phenomena. We argue that the quantum man- 
ifestations of both dynamical instability and equiparti- 
tion can be understood in terms of the depletion of the 
Floquet-Bogoliubov states, defined as the "low-energy" 
eigenstates of the evolution operator over one Bloch pe- 
riod. 

We furthermore go beyond the traditional single tra- 
jectory mean-field treatment and, following a recent sug- 
gestion [24j j . average the dynamics over an ensemble of 
trajectories given by the Husimi distribution of the ini- 
tial many-particle state. It is shown that this method 
is capable of describing important features of the many- 
particle dynamics. 

The paper is organized as follows: In Sec. [TT] we dis- 
cuss the mean-field dynamics, in particular the stabil- 
ity properties of the Bloch oscillation and its relation to 
chaotic dynamics. The corresponding many-particle sys- 
tem is analyzed in Sec. IIII1 mainly based on the Floquet- 
Bogoliubov states whose depletion properties provide a 
measure for the many-particle stability which can be 
compared to the mean-field behavior. We summarize our 
results and end with a short outlook in Sec. IIVI 



II. MEAN-FIELD DYNAMICS 

Evaluating the canonical equations of motion ([3]) with 
the Hamiltonian function ([2]) yields the mean-field equa- 
tions of motion of a BEC in a tilted optical lattice, the 
DNLSE: 

i hai = eim- ^(ai + i+ai-i) + g\ai\ 2 ai , ei=dFl. (4) 

Here ai(t) are the complex amplitudes of a mini BEC 
associated with the I th well of the optical potential, J 
is the hopping or tunneling matrix clement, d the lat- 
tice period, F the magnitude of the static force, and g 
the nonlinear parameter, given by the product of the mi- 
croscopic interaction constant W and the filling factor n 
(the mean number of atoms per lattice site). It should be 



noted that Eq. ((4]) can also be derived as a tight-binding 
approximation for the discretized Gross-Pitaevskii equa- 
tion. To simplify the equations we shall set the lattice 
period d and the Planck constant h to unity in the fol- 
lowing. 

Throughout the paper we shall use the gauge transfor- 
mation 

a,(t) -» ex P H(.9 + Fl)t]ai{t) (5) 

that eliminates the static term in Eq. Q. Note that the 
inclusion of g in the transformation is optional and is 
done to facilitate the stability analysis below. The effect 
of the static force then appears as periodic driving of the 
system with the Bloch frequency F: 

iai = -~ (e- iFt a l+1 + e+ iFt ai ^) + g (h| 2 - l) a t . 

An advantage of the gauge transformation is that one 
can impose periodic boundary conditions, ao(t) = cti(t), 
where we restrict ourselves to odd values of L. Equa- 
tion ([6]) also appears as a canonical equation of motion 
generated by the Hamiltonian function 

H ^ = -^E( e ^+i a ' + c - c -)+fEi a 'i 2 (i a H 2 -2) . 
1 1 

In this work, we shall be concerned mainly with al- 
most uniform initial conditions ai(0) ~ 1 which corre- 
spond to a BEC in the zeroth quasimomentum mode in 
the many-particle description in Sec. IIIII Note that we 
do not normalize |a;| 2 to unity. Of course, this initial 
condition is an idealization of the real experimental sit- 
uation, where initially only a finite number of wells are 
occupied. Nevertheless, an analysis of this situation pro- 
vides useful estimates which, in fact, can also be applied 
to the case of non-uniform initial conditions 25]. 

It is convenient to switch to the Bloch-waves represen- 
tation 

L 

b k = L- l/2 Y^e*v{iKl)ai , fc = 0,±l,...,±(Z-l)/2, 

1=1 

where k = 2nk/L is the quasimomentum (— ir < k < tt). 
After this canonical change of variables Eq. ^ takes the 
form 

ibk = — Jcos(k — Ft)bk 

+f E h t bl 2 b k3 Si% +k2 _ k3 -gb k , (9) 

kl ,/C2,&3 

where S^l, is the Kronecker 5 modulo L. For strictly uni- 
form initial conditions az(0) = 1, Eq. ([9]) has the trivial 
solution 

b (t) = VI exp (ij sin(Ft)) , ^o(i) = , (10) 



i.e. a Bloch oscillation with period Tb = 2ir/F. How- 
ever, it is well known that for g ^ the solution (flT)]) can 
be unstable with respect to a weak perturbation. The 
stability analysis of the solution p0[) . resulting in the 
stability diagram in the parameter space of the system, 
was presented in Ref. [13, H3|- We extend this analy- 
sis below, mainly following Ref. [22], using methods of 
classical nonlinear dynamics. 



A. Stability analysis 

From a formal point of view, the BO (|10|l is a periodic 
trajectory in the multi-dimensional phase space spanned 
by the cij or the b k , respectively. Using the standard 
approach, we linearize Eq. ([9]) around this periodic tra- 
jectory which leads to pairs of coupled equations, 



ib +k = - J cos(k - Ft)b +k + — \b \ b +k 



ib-k — — Jcos(k + Ft)b-b 



{blb*_ k 

jblb\ h (11) 



for k 7^ 0, where the initial amplitudes b±k(0) are arbi- 
trarily small. Substituting bo(t) from Eq. (flT)]) and inte- 
grating Eq. pip in time over n Bloch periods, we get 



b+k{t n ) 
6*fc(*n) 



= A?b! + A£b 2 



(12) 



Here t n = Tsn and A1.2 and bi2 are the eigenvalues and 
eigenvectors of the stability matrix 
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-r(t) -1 



dt 



(13) 



where the hat over the exponential function denotes time 
ordering, and 



/(*) 



f- 2J h 
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cos k] sin(Ft) 



(14) 



Note that the determinant - and therefore the product of 
the eigenvalues - of the (symplectic) stability matrix is 
equal to one. The trajectory is stable when both eigen- 
values lie on the unit circle but becomes unstable when 
they merge on the real axis and go in and out of the 
unit circle. In what follows we shall characterize this 
instability by the increment v, given by the logarithm 
of the modulus of the maximal eigenvalue: v = ln|Ai|, 
I Ai| > I A2 1 - The increment of the dynamical instability 
is parameterized by the quasimomentum k and, hence, 
there are (L — l)/2 different increments i^ fc ) > 0. For 
a stable BO, all of them should vanish simultaneously. 
Further details concerning the stability analysis can be 
found in Appendix lAl 

As an example we show the increment of the dynam- 
ical instability for the system with only L = 3 sites in 
dependence of F/J and gj J in Fig. [T] It is seen in the 
figure that the parameter space of the system is divided 
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FIG. 1: Increment of the dynamical (modulation) instability 
v as a function of the static force magnitude F/J and the 
interaction g/J for L — 3 sites. The yellow region in the 
lower panel corresponds to = for all k, i.e. the Bloch 
oscillation is stable. A magnification of the lower left corner 
is depicted in the upper panel and reveals a web of additional 
stability regions. 



into two parts by a critical boundary. If the number of 
sites is increased, the instability regions grow and the 
boundary approaches approximately the curve 
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F < 2.9 J 
F > 2.9 J 



(15) 



(see [Hj])- Figure [2] shows the stability diagram for L = 
63 lattice sites together with the boundary curve (fT5|) . 

In the right region of the diagrams v = v(F, g) = 
holds and BOs are always stable. Physically that means 
that a sufficiently large static force ensures stability even 
in the presence of nonlinearity. In the left hand side 
of the diagrams BOs are typically unstable but a closer 
inspection reveals a web of stability regions for smaller 
values of the parameters as can be seen in the upper 
panel in Fig. [TJ It should be noted that these stability 
regions are a particular property of systems with a small 
number of lattice sites. If L is increased, the regions 
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FIG. 2: Same as Fig. [1] (lower panel), however for L = 63 
sites. Also shown is the boundary given by Eq. (|15[l as a 
black curve. 




v {k) {F), k = 0, 1,... for J = 1, g = 0.1 and L = 5 sites 
(upper panel) as well as L — 63 sites (lower panel). With 
increase of the lattice size the instability regions for different 
k overlap and cover the whole region left to the critical line. 



where > overlap (see Fig. [3]) and for any point in 
the left side of the diagram there is at least one strictly 
positive increment of the dynamical instability. 

B. Relation to chaos 

In this section we investigate the implications of the 
stability on the dynamics of a classical ensemble of tra- 
jectories. We will address the relation between dynam- 
ical instability, classical chaos and the so-called self- 
thermalization, that is, equipartition of the energy be- 
tween the different quasimomcntum modes. Related 
questions have been studied in [26j in the context of one 



of the standard systems of classical chaos, the Fermi- 
Pasta-Ulam system (27}- There it was pointed out that 
a positive increment v of the dynamical instability is a 
necessary but not sufficient condition for the onset of 
developed chaos in a chain of coupled nonlinear oscil- 
lators. Here the term 'developed chaos' characterizes a 
situation where a chaotic trajectory explores the whole 
energy shell. Developed chaos also implies equipartition 
of the energy between the eigenmodes of the chain, a 
phenomenon often referred to as thermalization. 

The consideration of an ensemble of classical trajec- 
tories as compared to a single trajectory is well suited 
for understanding the systems behavior in the present 
context for two main reasons. First it is a convenient 
method to capture the generic features of perturbed ini- 
tial conditions due to an out-averaging of the special be- 
havior of an individual trajectory slightly differing from 
the BO (|10| . The second reason is that the mean- field 
description is an approximation in the spirit of a classi- 
cal limit of the many-particle system which has to be re- 
garded as the more fundamental description. This many- 
particle system, however, cannot be associated with a 
point in the (classical) mean-field phase space, but is 
rather equipped with a finite width due to the uncer- 
tainty principle, where the width decreases with increas- 
ing particle number. Thus the natural counterpart of 
the many-particle system within the mean-field descrip- 
tion is a phase space distribution rather than a single 
point. This phase space distribution can be conveniently 
replaced with a finite ensemble of classical trajectories for 
practical purposes. Further details of this considerations 
can be found in [24| and Appendix [B] As an example 
we depict in Fig. U the amplitudes a/, I = 1, . . . , L of 
an initial classical ensemble for the parameters N = 15, 
L = 5 considered in Sec. IIIII for 100 realizations. The 
histograms show the corresponding probability distribu- 
tions for the populations \ai\ 2 and the phases. 
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FIG. 4: Classical ensemble representing the many-particle ini- 
tial state for N = 15, L — 5 and 100 realizations. The char- 
acteristic width of the distribution is proportional to n _1,/2 . 

A convenient quantity for the characterization of the 
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Bloch dynamics in the higher dimensional classical phase 
space is the classical momentum given by 

= k E a ^ ai e ~ lFt - c A = E n 2 sin ( K - F< )- 

I k 

(16) 

For the periodic solution (fit)]) , that is, the BO, the mo- 
mentum oscillates in a cosine manner between —1 and 
1 with the Bloch frequency F. For neighboring initial 
values the behavior may differ considerably, depending 
on the stability of the BO: If the BO is unstable we ex- 
pect an exponential growth of the initial deviation con- 
nected to classical chaos, leading to thermalization. On 
the other hand, the naive expectation in the stable region 
is that a small initial deviation leads to a small deviation 
in the overall trajectory and therefore averaging over an 
ensemble does not change the Bloch oscillation behavior 
in principle. However, we are going to argue that this 
is only true in a certain range of the parameter space 
and indeed one can distinguish three regimes of classical 
motion instead of only 'stable' or 'unstable' behavior. In 
fact, we find that for large values of the static force F the 
ensemble average leads to a breakdown of the BOs and 
even a thermalization in the absence of classical chaos. 

Let us start our discussion with the unstable regime, 
that is, small values of F. The exponential growth of a 
deviation from the uniform initial condition is evidently 
related to chaos. In particular, as for the Fermi-Pasta- 
Ulam system [27] , for the system (J5J) positive increments 
of the dynamical instability are a necessary condition for 
the onset of developed chaos. To illustrate the impact 
on the Bloch dynamics we show the time evolution for 
a single trajectory of an ensemble mimicking an N = 15 
particle system in Fig.[5j The first panel in the figure de- 
picts the momentum (II 6[) and the second panel shows the 
mode populations |6fc(£) 2 | of the quasimomentum modes 
k as a function of time measured in units of Xj = 2ir/ J. 

For a single run, the momentum p(t) and the mode 
populations |6fc(i)| 2 start to oscillate irregularly after a 
transient time t v ~ ln(e/z/), e ~ n^ 1 / 2 , required for the 
modes with k ^ to take non-negligible values. These 
erratic oscillations are smoothed by averaging over the 
ensemble (consisting of 1000 trajectories in the present 
example) as shown in Fig. [5] In the ensemble average 
one observes a damped BO of p(t) and an equipartition 
of the mode populations converging to the values of X/L, 
i.e., a thermalization. It is clearly seen in the figure that 
the rate of thermalization actually determines the decay 
rate of the BO. In Sec. IIIII we will demonstrate that this 
(classical) mean-field dynamics agrees remarkably well 
with the full quantum many-particle behavior in view of 
the small number of three particles per site. 

One gets further insight into the relation to chaos by 
studying the finite time Lyapunov exponent 

\{t) = ln\Sa(t)\/t , (17) 

where Sa.(t) evolves in the tangent space according to the 



linear equation 

i^ = M[a(t)]*a(t) (18) 

(see Appendix [XJ). As an example, Fig. [7] shows the be- 
havior of X(t) for ten different trajectories from an ensem- 
ble (IB2|) . It is seen that X(t) converges to some constant 
value, so that the mean Lyapunov exponent (A) is a well 
defined quantity. Since the unstable periodic trajectory 
analyzed in subsection III Al is a member of this ensem- 
ble, the maximal increment of the modulation instability 
v = maxfci/ fc ' provides a reliable estimate for (A). Still, 
because the mean Lyapunov exponent depends on the en- 
semble (i.e. on the value of the filling factor n — N/L), 
generally (A) =/= v. In particular, (A) is found to be a 
smooth function of F, while the maximal increment of 
the modulation instability is a non-analytic function of 
F with discontinuous first derivative. As a consequence, 
when we cross the critical line in the stability diagram, 
the system shows a smooth transition from the regime of 
decaying BOs to the regime of persistent BOs. (For ex- 
ample, for the parameters of Fig. [6] this change happens 
in the interval 0.1 < F < 0.4.) 

We now turn to the parameter regime of stable BOs, 
for larger values of F where the increment is zero. Here 
we will distinguish two different types of behavior in the 
ensemble average: Persistent BOs for intermediate val- 
ues of F and decaying BOs connected to thermalization 
in the limit of large F. The regime of persistent BO is 
shown in Fig. [8] for an example with F = 0.4. No energy 
exchange between quasimomentum modes is seen which 
means that the system dynamics is at least locally regu- 
lar. In other words, for the considered moderate F the 
periodic trajectory (|10p is surrounded by a stability is- 
land. Moreover, the size of this stability island should 
be large enough as compared to the characteristic width 




FIG. 5: Mean-field BO of the momentum (upper panel) and 
the dynamics of the population of quasimomentum modes 
(lower panel). These results for a single trajectory oscillate 
erratically. Parameters are L = 5, J = 1, g = 0.1, and 
F = 0.1. 
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FIG. 6: Mean-field BO of the momentum (upper panel) and 
dynamics of the populations of quasimomentum modes, aver- 
aged over an ensemble of 1000 trajectories of initial conditions 
(|B2|) . The other parameters are the same as in Fig. [S] An 
equipartition between the quasimomentum modes results in 
the decay of BO. 

of the distribution (|B2[) . so that the majority of the tra- 
jectories are stable. In this case, the ensemble averaging 
is of little influence, resulting only in a small decrease of 
the amplitude. 

The behavior for large values of F and its relation 
to classical chaos is more surprising: An analysis of the 
phase-space structure of the system (O reveals the vol- 
ume of the regular component to grow with F and for 
F — > oo the system becomes integrable [23|. (This inte- 
grable regime has already been observed in [281].) Thus, 
one would first expect persisting BOs. However, the ob- 
served behavior is quite different. As an illustration we 
show an example for an individual trajectory in Fig. [SJ 
Here one observes a quasiperiodic behavior. When aver- 
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FIG. 7: Finite time Lyapunov exponent X(t) for ten different 
trajectories from the ensemble (|B2[I . (Same parameters as in 
Fig.®) 



aged over an ensemble this leads to a decay of the BO as 
depicted in Fig. [TU1 This behavior can indeed be under- 
stood analytically. In the limit F — > oo the populations 
of the lattice sites are frozen and the amplitudes ai evolve 
only in the phase according to 

ai(t) woi(0)exp[-* 5 (|oi(0)| a -l)t] • (19) 

The evolution (p~9j) for the amplitudes ai(t) immediately 
implies the observed quasiperiodic dynamics of the am- 
plitudes bk(t) (see Fig. [9J. It can be shown that the 
decay due to the dephasing arising for the ensemble av- 
eraged dynamics obeys an exp(— r y r t 2 )-\&w, where the 
coefficient 7 r is proportional to g and inversely pro- 
portional to n. Remarkably, the quasiperiodic dynam- 
ics also implies an equipartition between the quasimo- 
mentum modes (see Fig. [9]) in the absence of classical 
chaos. However, there are characteristic differences com- 
pared to the decay and the thermalization processes in- 
troduced by classical chaos: First, the relaxation con- 
stant 7 r for the dephasing decay depends on the charac- 
teristic width of the distribution function (|B2j) and de- 
creases as 1/n in the semiclassical limit. On the con- 
trary, the relaxation constant 7 C for the chaotic decay 
decreases only as 1/lnfi. Second, the chaotic decay fol- 
lows (p(t)} — exp(— 7 C £) sin(wsi), while for the dephasing 
decay we have (p(t)} = exp(—j r t 2 ) sin^si). 

Going ahead we note that both the regimes of sta- 
ble BO and the one of decaying BO in the presence of 
classical chaos closely resemble the full many-particle dy- 
namics when averaged over an ensemble. However, in 
what follows we shall find that there is an additional 
many-particle feature in the dynamics in the third regime 
of large F. Here the decay of the Bloch oscillations is 
present in the many-particle system as well and for short 
times the mean-field ensemble even quantitatively resem- 
bles the many-particle dynamics. However, the many- 
particle system shows a periodic revival behavior as al- 
ready pointed out in [28| which, as a pure quantum phc- 
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FIG. 8: Same as Fig. [6] however for F = 0.4 showing a stable 
Bloch oscillation. 
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FIG. 9: Same as Fig. [5] yet for F = 10 in the strong field 
region. 



nomenon, cannot be captured by the mean-field ensem- 
ble. 



III. MANY-PARTICLE DYNAMICS 

In this section we study the many-particle counter- 
part of the mean-field system discussed in the preceding 
section. This is described by the driven Bose-Hubbard 
Hamiltonian 

= -f E (f'ti+i* + h - c ) +yI>(fti-2). 
i i 

(20) 

We also confine the system to L lattice sites (L chosen to 
be odd) and apply periodic boundary conditions. Then 
the Hilbert space of system (|20|) is spanned by the Fock 
states |n) = \rii, ■ ■ . , iil), where J2i ni = N is the to- 
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FIG. 10: Same as Fig. [6] yet for F = 10 in the strong field 
region. 



tal number of atoms. The dimension of the Hilbert space 
is equal to ^f7^rw • Using the canonical transformation 

b k = L~ x l 2 J^i exp(i27rfcZ/L) a./, the Hamiltonian (|20[) can 
be presented in the form 

H = -J^cos(k- Ft)b\b k 

h 

+ ^ 2^ b ki b k2 b k3 b ki 4i+fe 3 ,fe 2 +fe4 ) ( 21 ) 

where we omit the constant term W^ k b k b k - The 
basis vectors of the Hilbert space are now given by 
the Fock states in the quasimomcntum representation: 

\n-(L-i)/2,---,n^i,no,n + i,...,n (L _ 1 y 2 )- In the co- 
ordinate representation the Fock and the quasimomen- 
tum Fock states are given by the symmetrized product 
of the Wannier and Bloch functions, respectively. Here 
we are interested in the solution of the time-dependent 
Schrodinger equation with the Hamiltonian (|2~lj) for ini- 
tial conditions given by a BEC of atoms in the zero quasi- 
momentum state, i.e. |$g) = I • • ■ i 0, N, 0, . . .) q . This is, 
in fact, equivalent to an SU(L) coherent state, namely 

'* 0> = vM^lj'^ = ^ Cn|n> (22) 

where |n) with n = (ni,...,rii) is a Fock state and 
c n = ^ji^'.^ j. In general the SU(L) coherent states 

are equivalent to the fully condensed states, where our 
special choice approximately corresponds to the ground 
state of the system at F = if the condition W <C J is 
satisfied. 



A. Floquet-Bogoliubov states 

First we address the question of a manifestation of the 
dynamical instability in the many-particle quantum sys- 
tem. It is argued below that, similar to the static case 
F = 0, the quantum counterpart of the dynamical in- 
stability is Bogoliubov's depletion of the condensate . 
We begin with an alternative derivation of the common 
Bogoliubov spectrum for F = [3(| which we shall then 
adopt to the case F ^ 0. 

For an infinite number of particles, the Bogoliubov 
states can be constructed from |*P ) by applying the de- 
pletion operators 

oo 

^^^f^lM)". (23) 

n=0 

which transfer particles from the zero quasimomentum 
state to the states ±fc : 

|*> = n 5(fc) l*o) ■ (24) 

k>0 
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In Eq. ((23)) the coefficients c^f* should be determined 
self-consistently so that the wave function |^) satisfies 
the stationary Schrodinger equation with the Hamilto- 
nian (ED with F = 0, 



H(F = 0)\V) =E\V) . (25) 
, (f2"4"f are equivalent to the ansatz 



Equations 



\^ k) ) = J2c^\..,n,..,N-2n,..,n,..) q 



(26) 



n=0 



where n particles are redistributed from the zero quasi- 
momentum state to the state k and to the state —k. Note 
that here, as well as in ([23]) , (|24|) . the limit N — ► oo, 
W — > with constant g — WN/L is assumed, which jus- 
tifies a factorization of the eigenvalue problem (f25l) into 
(L — l)/2 independent eigenvalue problems. Substitut- 
ing the ansatz (|26|) into (|25|) and taking the above limit 
the original problem reduces to the diagonalization of a 
tridiagonal hermitian matrix A^ k ' with matrix elements 



m=2(g + 8)n 



A 



g(n+l) 



(27) 



where 5 = J (1 — cos n). The spectrum of the matrix A^ k ' 
is equidistant with a level spacing given by the Bogoli- 
ubov frequency — 2y2g6 + 5 2 . Correspondingly, 

the eigenvectors c( fc ) of A^ define the Bogoliubov states 
of the Bose-Hubbard model. It is worthwhile emphasiz- 
ing that for a finite N these Bogoliubov states provide 
only an approximation to the actual eigenstates of the 
Bose-Hubbard system. While this approximation can be 
quite accurate for the ground and low-energy states, it 
fails for the high-energy states [30|. The BEC depletion 
of these states, defined by 



Nr 



E 

k>0 n=l 



5>i 



,W|2 



(28) 



provides a quantitative criterion for the validity of the 
Bogoliubov approach. Namely, if No << N then the 
approach is justified. On the contrary, No ~ N (which 
is the case for high-energy states) means that the actual 
structure of the eigenstates has nothing to do with the 
presumed Bogoliubov structure (|24p . Furthermore, it is 
convenient to order the states according to the depletion 



Now we are in the position to discuss the case F ^ 0. 
Since we are interested in the Bloch dynamics, it is useful 
to introduce the evolution operator over one Bloch period 

fil, 



U = exp 



H{t)dt 



(29) 



where H (t) is given in (|2"Tj) and the hat over the exponen- 
tial function again denotes time ordering. We are looking 
for those eigenstates of U, 
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FIG. 11: Left panel: Coefficients of the first three Bogoli- 
ubov states I^P} = Cn\n, N — 2n, n) of the three-sites Bose- 
Hubbard model. Right panel: Floquet-Bogoliubov states. 
(Parameters J = 1, g = 0.1 and F = 0.4.) 



which have the Bogoliubov structure Substituting 
the ansatz ([26]) into ([30]) we obtain an equation for the 



coefficients c^f 1 : 



exp 



-i [ B A {k) {t)dt 
Jo 



,(*0 



.•2^\„( k ) 



exp (-*— j C W , (31) 



U |*) = exp{-i2TrE/F) |*) , 



(30) 



where the diagonal elements of the matrix A^ k \t) are 
now given by 

Ai k Ut) = 2[g + 6cos(Ft)]n, 6 = J(l - cos k) . (32) 



As an illustration the right panel in Fig. [TT] shows the 
coefficients of the first three Floquet-Bogoliubov states in 
the quasimomentum Fock basis for L = 3, g = 0.1, and 
F = 0.4 [U, where the depletion is N D = 0.113, 2.341, 
and 4.568 particles, respectively. For the sake of compar- 
ison, the left panel in the figure shows the Bogoliubov 
states for the same value of g = 0.1, where the depletion 
is N D = 0.002, 2.006, and 4.010. It is seen in Fig. [TT] 
that the depletion of the driven BEC is larger than the 
depletion of a stationary BEC which was found to be a 
typical situation in further numerical studies. 

The depletion l|28[) provides useful information on the 
BEC stability. For the 'ground' Floquet-Bogoliubov state 
the dependence No — Nd(F) is depicted in Figs. [T2l and 
[T~3l It is seen in Fig. [12] that, as a function of F, No 
diverges at the points where the classical increment of 
the dynamical instability v (also shown in the figure) 
takes positive values. (More precisely, the depletion can- 
not be larger than N. However, since the Bogoliubov 
theory refers to N — oo, it can be formally considered 
as infinite.) This 'divergence' of No means that in the 
regions of instability the eigenfunctions of the evolution 
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FIG. 12: Number of the depleted particles Nr> (red curve) as a 
function of the field strength F (logarithmic scale) compared 
to the classical increment of dynamical instability v — v(F) 
(blue curve). Parameters are J — 1, g — 0.1, L — 3. The 
upper panel magnifies the region of small F in the lower panel. 

operator differ considerably from the Bogoliubov struc- 
ture. In fact, they are chaotic in the sense of quantum 
chaos [13| . It is also seen in the figure that the deple- 
tion increases linearly for F — > oo. Thus the discussed 
Floquct-Bogoliubov states cannot be eigenstates of the 
evolution operator in the limit of large F, although the 
classical increments v^> vanish identically for F — ► oo. 
We come back to this point in the next subsection. 

The intervals of small depletion observed in Fig. [T2l 
depend, of course, also on the interaction g. The full 
parameter dependence of the depletion (|28[) is shown in 
Fig. [13] For relatively low values of gj J and Fj J (upper 
panel) we find a highly organized web of stability regions, 
whereas the behavior for larger parameters (lower panel) 
shows a simpler structure. These many-particle results 
can be directly compared with the mean-field stability di- 
agrams in Fig. [1] which confirms the relationship between 
BEC depletion and mean-field stability. 

For the sake of completeness we also briefly discuss 
the corresponding quasienergies E (see Eq. ([50])) which 
are defined modulo h times the Bloch frequency, i.e. 
E a ,j = E afi + jF, j = 0, ±1, .... The lower panel of 
Fig. [14] displays the quasienergies of the five Floquet- 
Bogoliubov states with smallest Nd . The system param- 
eters are J = 1, g = 0.1 and L = 3 as in Fig. [T2J For 
clarity, three Floquet zones are shown in the figure. As 
expected, the spectrum is equidistant with a level spac- 
ing correlated to the depletion. Note the apparent irreg- 
ularity in the windows around F = 0.2 and F = 0.17, 
which is related to the region of dynamical instability of 
the mean-field system. For comparison, the upper panel 
shows the number of depleted particles No as in Fig. [T^] 
with a linearly scaled _F-axis. In the regions of finite de- 
pletion the most stable quasienergy states appear to be 
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FIG. 13: Number of the depleted particles No as a function 
of the static force magnitude Ff J and the interaction g/ J for 
L = 3 sites. The depletion is very small in the yellow region 
of the lower panel. The upper panel magnifies the lower left 
corner of the lower panel and reveals additional regions of 
small depletion. 



very sensitive against a variation of F. 



B. Bloch oscillations 

The microscopic dynamics of BOs was considered ear- 
lier in a number of papers summarized in Ref. [1 31 ] with 
focus on the regime of low filling factors n = N/ L < 1 . 
In the present work, to make a link to the mean-field 
dynamics, we simulate BOs for a relatively large filling 
factor. 

We investigate the dynamics in dependence on the 
force F for N — 15 and N = 20 atoms in a lattice with 
L = 5 sites. Since the plots are very similar, only the 
results with iV = 15 are shown. The value of the micro- 
scopic interaction constant is set to W — 0.1/n, so that 
the macroscopic interaction constant is g = 0.1, and the 
hopping matrix element ,7=1. The initial wave function 
is chosen as the SU(L) coherent state (|2"2"j). 
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FIG. 14: Lower panel: Quasienergy spectrum of the Floquet- 
Bogoliubov states as a function of the field strength F (same 
parameters as in Fig. I12[) . Upper panel: Number of the de- 
pleted particles Nd- 



FIG. 15: Bloch oscillations of N — 15 atoms in a lattice with 
L = 5 sites for J = 1, W = 0.1/3, and F = 0.1 (top), F = 0.4 
(middle), and F — 10 (bottom). Shown is the many-particle 
mean momentum p(t) given in (|33|) . 



In Fig. [T5] we show the dynamics of the many-particle 
counterpart of the mean- field quantity p^|) . that is, the 
expectation value of the many-particle momentum 



Pit) 



1 

2iN 



ai e 



-iFt 



h.c. 



m) (33) 



for N — 15 and the three values of F chosen in different 
dynamical regimes as shown in Figs. l6l [8l and [T0l 

The upper panel of Fig. [TBI corresponds to F = 0.1, 
which falls into the region of dynamical instability of the 
mean-field system. Here the quantum many-particle BO 
decays in very good agreement with the mean-field en- 
semble shown in the upper panel of Fig. [5] (note that finer 
details are also reproduced) . This demonstrates that the 
ensemble averaged mean-field dynamics is capable of de- 
scribing important aspects of the full many-particle sys- 
tem [241 ]. 

In the middle panel of Fig. [15] we have F — 0.4, where the 
system is stable. Here the quantum many-particle BOs 
persist in time, also agreeing with the ensemble-averaged 
mean-field dynamics. As discussed above, in this regime 
the ensemble average is of little influence, i.e. the state is 
fully condensed and can be described by a single mean- 
field trajectory instead. 

The third case, F = 10, depicted in the lower panel 
of Fig. [T5l requires a separate consideration. Indeed, as 
mentioned in the previous subsection, in the limit of large 
F the Floquet-Bogoliubov states are not eigenstates of 
the evolution operator (|2l)|) . Instead it can be shown 
that those are the Fock states |n) [l3j]. Thus the time 
evolution of the wave function is given by 



1*0)) = 5] Cnex P 



i=i 



(34) 



Equation (|34| implies a periodic recovering of the initial 
state (f2"2"| at times which are multiples of Tw = 2tt/W 
and, hence, periodic revivals of BOs 28]. It should be 
stressed that these revivals are a pure quantum many- 
particle effect due to the finiteness of W and n. This 
constructive interference cannot be explained within the 
ensemble-averaged mean-field approach. However, the 
breakdown can be described by the ensemble averaging 
and is due to dephasing, as explained in Sec. Ill Bl 



IV. SUMMARY 

We have studied an TV-particle system, a Bose- 
Hubbard Hamiltonian with linearly increasing on-site en- 
ergies. This system can be conveniently reduced to a fi- 
nite lattice with L sites by using gauge transformation 
and imposing periodic boundary conditions. 

Such a model can be used to describe important fea- 
tures of realistic systems, as for instance the dynamics 
of cold atoms or BECs in an optical lattice under the 
influence of the gravitational field [H, i, i, i, i, @, 0, 0, i, 
ITol . [Til . [l2| , or many-particle systems in ring shaped op- 
tical lattices as proposed in [321 ] with additional driving. 
It should be noted, however, that it neglects a number 
of features. The space dimension is reduced to a quasi 
one-dimensional setting, a decay of the system via Zener 
transitions to higher bands is excluded, the lattice is dis- 
cretized and truncated, and the interaction is simplified. 
Nevertheless, this model has been found to describe ex- 
perimental results quite well. Further, the Bose-Hubbard 
model is of interest in its own right, as evident from the 
large number of studies exploring its properties which are 
remarkably rich. 

In the present paper we have studied the Bose- 
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Hubbard model for relatively small number of lattice sites 
L < 5 but relatively large number of particles up to 
TV = 20 to make a link with the mean-field dynamics. 
Our aim was twofold: First, we have demonstrated how 
the modulation instability observed in the mean-field sys- 
tem are manifested in the many-particle case. We have 
shown that a reasonable measure of the A^-particle insta- 
bility is provided by the quantity (|28p . Second, we have 
explored the possibility to describe the time-evolution of 
many-particle expectation values in terms of mean-field 
trajectories, averaged over an ensemble constructed from 
the ST/(L)-phase space distribution of the initial many- 
particle state. We found that this (averaged) mean-field 
dynamics agrees remarkably well with the full quantum 
many-particle behavior in a number of cases. The only 
exception was the presence of many-particle revivals for 
strong fields which is a pure quantum phenomenon. 

These observations suggest an application of the mean- 
field ensemble method to investigate the properties for 
larger lattices and larger particle numbers where many- 
particle computations are much more difficult or virtu- 
ally impossible. Furthermore, it will be of interest to 
study the interrelation between classical chaotic motion 
and stable or decaying Bloch oscillations for more lattice 
sites where one can possibly make contact with recent 
related studies of the force free, F = 0, case, both for 
the mean- field [33l l3~i| and the many-particle description 
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where M [a(i)] is a 2L x 2L matrix of the following struc- 
ture: 



M[a(t)] 



A + gB gC 
-gC* -{A + gB)* 



(A3) 



Ai <m = - ^ (Si+i,me lFt + 8t-i <m e lFt ) , (A4) 
Bi, m = \ai{t)\ 2 5u m , C^ m = Clf(t)5i, m . (A5) 

Inserting the trajectory (IA1|1 into (|A2[) , the linear equa- 
tion (|A2|) takes a form where the matrix M(t) is periodic 
in time. Finally, we introduce the stroboscopic map for 
the discrete time t n — Tqii: 



6a(t n+1 ) = USk(t n ) . U = exp 



-I j M{t)dtj 



(A6) 

This stability matrix U is symplectic and hence the con- 
sidered periodic trajectory (|Aip is stable if and only if all 
its eigenvalues lie on the unit circle. 

Using the unitary transformation U — > VUV" 1 , where 



V = 



T 
T 



rp _ r-1/2 
J-l.rn — L 



exp 



i2itk 



(A7) 



the stability matrix U can be factorized into L decou- 
pled 2x2 matrices which can be labeled as with 
k = 0,±1,...,±(L- l)/2. The matrix U^> is not of 
interest here because its eigenvalues are always located 
on the unit circle. The explicit form of the matrix U^ ±k ^ 
is given by Eq. (fT3j) and, depending on F, its eigenvalues 
Xi.2 either lie on the unit circle (with A2 = A*) or on 
the real axis (with A2 = 1/Ai). We have stability in the 
first, and instability in the second case. 
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APPENDIX A: STABILITY MATRIX 

It is instructive to consider the DNLSE © from the 
viewpoint of the general theory of nonlinear dynamics. 
Then the solution with the initial condition a/(t = 0) = 1 , 

ai(<) =exp(i^sin(.Fi)) (Al) 

is nothing other than a periodic trajectory in a 2L- 
dimensional phase space. Thus one may address the 
question of stability of this periodic trajectory. 

Denoting by da = (Sai, . . . , Scil, 6a*, . . . , 6a* L ) T the de- 
viation from an arbitrary trajectory &(t) and linearizing 
the DNLSE around this trajectory we have 

i ^6a = M [a(t)] 6a , (A2) 



APPENDIX B: MEAN-FIELD ENSEMBLE 

Let us recall that the mean-field dynamics of the site 
amplitudes a; (or, alternatively the quasimomentum am- 
plitudes bk) appears as a canonical Hamiltonian evolution 
in a 2L-dimensional complex phase space which conserves 
the norm, i.e. the phase space is the surface of a com- 
plex sphere. In order to approximate the many-particle 
dynamics, we construct an ensemble of mean- field tra- 
jectories representing the initial many-body state |^(0)). 
This is achieved most conveniently by using the quantum 
Husimi phase space distribution [24| . the projection onto 
SU(L) coherent states which in the quasimomentum rep- 
resentation is given by 




with b = (6_(i_i)/ 2> -.,&o,»,6(L-i)/2) and EfcN 2 = 
L. For the initial many-particle state (f2"2"]) this Husimi 
density Q(h) = |(b|5'(0))| 2 is easily calculated as 

Q(b)=B\b \ 2N , (B2) 
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where B is a normalization constant. Note that the prob- 
ability density for the vector components depends only 
on the zero-mode probability |&o| 2 and is strongly local- 
ized in the region close to its maximum value |6o| 2 = L 
for large particle number N . 

Numerically, one can construct the ensemble (|B2|) by 
means of a rejection method [36|: First one generates L 
randomly distributed real numbers in the unit inter- 
val with random phases 4>k and normalizes the random 



vector with components — e l ^ k to unity. Another 
random number v in the unit interval is chosen in order 
to decide if the generated vector is accepted as a part of 
the ensemble: it is accepted if v < \zq\ 2N and rejected 
otherwise. Finally we renormalize as bk = y/Lzk- If de- 
sired, the ensemble can be Fourier transformed to yield 
the corresponding ensemble of lattice site amplitudes ai . 
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